fig.path <- "results/graph/"
fig.name <- "pap_fig_03.pdf"

yvars <- c("par.onl.sh")

dist <- c()
effect <- c()
ses <- c()

steps <- seq(500, 10000, 500)
preserve <- Re
for (i in steps){
  j = as.character(i/100)
  k <- str_pad(j,2,pad = "0")
  Tr = paste0("march",j)
  fTr = paste0("factor(", Tr, ")")
  print(Tr)

  ftreat <- paste0("factor(treat_", k, ")")
  treat <- paste0("treat_", k)
  fcontrols <- paste0("dpost:",controls)
  vars <- c( paste0("factor(dpost)*", ftreat),"factor(dpost)" ,ftreat)
  outvars <- c( paste0("factor.dpost.TRUE.factor.",treat,".TRUE"),
    "^factor.dpost.TRUE$" ,
    paste0("^factor.", treat, ".TRUE$"))
  labvars <- c("DPost X March","DPost", "March")


  m1 <-felm(as.formula(RegFor( y = yvars[1] , x = c(vars, fcontrols) ,
      FE = "county" , IV="0", clust = "klust" )),
           data= Re)

  beta <- coefficients(m1)[paste0("factor(dpost)TRUE:", "factor(", treat , ")TRUE")]
  se <- m1$cse[paste0("factor(dpost)TRUE:", "factor(", treat , ")TRUE")]
  dist <- c(dist, j)
  effect <- c(effect, beta)
  ses <- c(ses, se)
}
kol <- ghibli_palette(name = "PonyoMedium")[c(4)]
Re <- preserve
coefs <- data.frame(dist, effect, ses)
coefs$up10 <- coefs$effect + 1.645*coefs$se
coefs$lo10 <- coefs$effect - 1.645*coefs$se
coefs$up5 <- coefs$effect + 1.9599*coefs$se
coefs$lo5 <- coefs$effect - 1.9599*coefs$se
coefs$dist <- as.numeric(coefs$dist)/10

p.t <- ggplot(data = coefs, aes(x = dist)) +
geom_pointrange(aes( y= effect, ymin=lo5, ymax = up5) , 
  size = 1, fatten = 2,  alpha = 0.9, 
  position = position_dodge(width = 0.4),
  col = kol) +
geom_linerange(aes(ymin=lo10, ymax = up10) , 
  size = 3,  alpha = 0.35, 
  col = kol,
  position = position_dodge(width = 0.4)) +
geom_vline(xintercept=0,linetype="dotted", size=0.2) +
geom_hline(yintercept=0,linetype="dotted", size=0.2) +
theme_light(base_size = 16) + labs( x = "Distance parish - march (km)" , y = "Effect DPost x March") +
guides(shape = guide_legend(override.aes = list(size = 0.5)))

ggsave(plot = p.t , paste0(fig.path, fig.name), width = 8, height = 4.5)
rm(list=
    setdiff(ls()[sapply(ls(), function(x) any(class(get(x)) == 'data.frame'))],
    c("Re","Ind")))